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ABSTRACT 

Given the importance of simulating hydromagnetic processes that impact star 
formation, we have earlier developed a 3D adaptive mesh approach that allows us to 
include hydromagnetic processes during the formation and evolution of cores, discs, 
and stars in observed regions of star formation. In this paper, we take the next step 
in this program - namely - to d evelop a modified version of the 3D adaptive mesh 
refinement (AMR) code FLASH (|Frvxell et al.ll2000l ) in which the ambipolar diffusion 
of the magnetic field in poorly ionized molecular gas is implemented. We approach 
the problem using a single-fluid approximation to simplify numerical calculations. In 
this paper, we present a series of test cases including oblique isothermal and non- 
isothermal C-shocks. We also present a study of the quasi-static collapse of an initial 
uniform, self-gravitating, magnetized sphere that is initially supported by its magnetic 
field against collapse (i.e. magnetically subcritical). Applications to the collapse of a 
pre-stellar Bonnor-Ebert sphere are presented in a companion paper. 

Key words: MHD - Shockwaves - methods: numerical - stars: formation - ISM: 
clouds - ISM: magnetic fields. 



1 INTRODUCTION 

A powerful new theoretical picture of star formation, 
known as turbulent fragmentation, has emerged over the 
last decade. It emphasizes the role of supersonic turbu- 
lence in producing the filamentary structure and mass 
spectrum of overdense regions in Giant Molecular Clouds 
(GMCs) that ultimately collapse to form clusters of stars. 
A large number of simulations from a number of differ- 
ent groups have played a decis ive role in the evolution 
of these ideas (e.g. review s by lElmegreen fc Scald |2004 
Mac Low fe Klessenl |2004|; I Allen et alj 120071 : iKlein et all 
20071 . iBonnell. Larson, fc Zinneckerl l2007r i. It has been 
demonstrated that 3D turbulence provides a reasonably 
good explanation of the origin of the core mass function 
(CMF). Cores are the overd ense regions in which th e 'mi- 
crophysics' of star formation (|McKee fc Qstr ikcr 2007) - the 
formation of single or binary stars by gravitational collapse - 
plays itself out. F its of numerical mod els to observed CMFs 
are encouraging IjPadoan et al.l 11999 ). Another important 
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aspect of this picture is that the vorticity that is gener- 
ated in such oblique shocks naturally explains the origin 
of angular momentum and format ion of accretion discs ob- 
served around all young stars (e.g. Ijappsen fc Klessenll2004l : 
iTillev fc Pudritj 120041 ). The subsequent fragmentation of 
discs depends both on their mass and angular momentum 
and is import ant for the genesis of binary stellar systems 
(|Toomrdll964f ). Only with the advent of fully 3D simulations 
has it been possible to move beyond highly idealized spheri- 
cal or axisymmetric modeling to investigate non-equilibrium 
structures that bear close resemblance to those that are ob- 
served. 

Magnetic fields can significantly affect these processes 
on both the macroscopic and microscopic scales of star 
formation. Although only a few 3D numerical hydromag- 
netic simulations are available, the results show that if 
molecular clouds have magnetic energy densities close to 
gravitational self-energy, fragmentation into stars can be 
strongly affected. The so-called mass-to-flux ratio, F = 
M/$ = 27rG 1//2 E/B, controls the evolution of clouds in 
in the limit that their magnetic fields are strongly cou- 
pled to the gas ('ideal magnetohydrodynamics (MHD)'). 
For values of V ~ 1 — 3, magnetic fields force the col- 
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lapse of structure into pancakes and can strongly limi t 
the fragmentation into cores (e.g. iTillev fc Pudrit3l2007l) . 
Robust fragmentation of overdense regions within GMCs 
into CMFs requires th at they are sufficiently magnetically 
supercritical (r > 1) jHeitsch. Mac Low, fc Klessenl l200ll ; 
lOstriker. Stone, fc Gammiell200lh ~ It is important to real- 
ize, however, that the turbulent fragmentation of supercrit- 
ical clouds still produces cores that have a distribution of 
magnetizations, most of which are more flux loaded (lower 
mass to flux ratio) than their large scale surroundings. In 
fact, even a strongly supercritical clo ud can produce some 
cores that are critically ma gnetized ()Padoan fc Nordlundl 
ll999l ; lTillev fc Pudritzll2007h . This agrees with the observed 
presence of dynamically significant m agnetic fields in a sig- 
nificant number of cores (but not all: ICrutcherll 19991 ). 

The strength of these magnetic effects can be much re- 
duced, in principle, if fields are initially very weak or they 
are poorly coupled to the denser gas. The high column den- 
sity of molecular gas implies that the gas-fiel d coupling is 
far from perfect l|Umebavashi fc N akano 198l]). In the bulk 
of the gas in molecular clouds, we observe ionization abun- 
dances on the order of (10 _7 )riH. At the average densities in 
cluster forming gas within molecular clouds (e.g. densities of 
~ 10 4 cm" 3 ), the ionization fraction of molecular gas is lower 
but still leaves the field somewhat coupled and active in the 
dynamics. Only the charged species in a gas are directly cou- 
pled to the field, and it is through elastic collisions with neu- 
tral particles that the gas as a whole responds to magnetic 
forces. Provided ion-neutr al collision rates are large enough 
l|Li. McKee. fc KleinllioOrl l. charged species can attain sig- 
nificant velocities differing from the neutral gas velocity as 
they are dragged by the magnetic field lines in a process 
called ambipolar diffusion (AD). This process is also known 
as ambipolar drift, or as ion-neutral drift by plasma physi- 
cists. With increasing slippage, this process can generate 
significant amou nts of thermal energy - more so tha n cos- 
mic ray heating (jPadoan. Zweibel. fc Nordlu nd 2000) - and 
can act as a significant mechanism for the damping of tur- 
bulence in molecular clouds. In order to understand some of 
the effects of this complicated process, it is helpful to think 
of AD as lying somewhere in between the purely hydrody- 
namic (HD) behaviour of gas, and the perfectly coupled limit 
in which magnetic forces can strongly reshape purely HD 
behaviour. 

The modeling of AD has been an i mportant ingredi- 
ent i n star formation theory for de cades ([Mestel fc Spitzerl 
Il95fj ; IShu. Adams, fc Lizandfl987l ). This process was first 
used in an astrophysical context bv lMestel fc Spitzerl ()l956l ) 
to describe the effective loss of magnetic support in a mag- 
netically subcritical (r < 1) cloud of gas that is other- 
wise thermally supercritical, leading to its eventual col- 
lapse. Many early theoretical treatments of AD in the liter- 
ature involve idealized one or two dimensional treatments in 
which only the evolution of non-rotating spheres or cylin- 
ders are treated. The role of ambipolar diffusion has been 
studied in the context of qu asi-static collapse of magnet- 
ically supported thin disc s (Fiedler & Mousch oviasl 1 19931 ; 
iBasu fc Mouschoviasl Il994l ; ICiolek fc Basul 120061 ). allowing 
for simplified geometries. This has recently been advanced 
by the correlated evolution of chemical species in the collapse 
dDesch fc Mouschovi as 2001; Naka no. Nishi. fc Umebavashil 
120021 ; iTassis fc Mouschoviasll2007h . ~ 



Simulations using the 3D Zeus-MP static grid code 
have studied the structure and stability of magnetic 
shocks in the presence of strong ambipolar diffusion (C- 
shock s ; lTothlll994lMac Low et a\lll995l; ISmith fc Mac Low! 
IMac Low fc Smithl Il997l ; \Li et all I2006T ). includ- 
with ambipolar diffusion 
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|Mac Low fc Smithll 19971 ). IMac Low et all (| 19951 ) started by 
implementing a single-fluid approximation, later develop- 
ing a t wo-fluid model in order to capture shock insta- 
bilities l|Mac Low fc Smithl Il997l ). Single-fluid approaches 
have been use d to study star-formation and astrophysi- 
cal turbulence ( Hosking & Whitworth 2004b; Padoa n et al. 
l2000l;IZweibel fc Brandenburdl 19971; iFiedler fc Mouschovias 



ll993j). Multi-fluid codes have recently been developed to 
study the effects of ambipolar diffusion on turbulence 
and the formation of star forming clumps and shocks 
|Oishi fc Mac L ow 2006; X/Tet al.ll2006l). including more gen- 
eraliz ed physics (|Falld [20031 ; lO'Sullivan fc Downed 120071 . 
2006) all of which implement methods that avoid small 
timesteps. 

In this and a series of upcoming papers, we investi- 
gate critical processes in star formation taking the effects 
of ambipolar diffusion into account through the use of the 
single-fluid approximation (carried out in £|2.3|1 . Given cer- 
tain physical conditions present in star-forming regions, we 
can ignore the charged fluid's inertia and pressure in a gas. 
This allows one to then treat only the dynamics of a single 
neutral fluid, with additional effects of the magnetic field 
due to elastic collisions with the charged fluid. This assump- 
tion is typically made for the simplicity in implementing the 
approach and in the numerical calculations. We emphasize 
that this approximation is physically valid in star forming 
regions where we have low ionization and strong coupling of 
ions with the field. However, it must be noted that in small 
regions in molecular clouds this method is expected to be 
limited. As an example, the small central regions within col- 
lapsing cores have densities exceeding n n > 10 10 cm -3 such 
that grains will play a dynamical role. This is also true in 
C-shocks (shocks mediated by ambipolar diffusion), where 
the width of the shock will be affected by the presence of 
grains. The single-fluid approximation and its limitations 
are further discussed in £12.21 where we show that it is still 
a useful approximation. 

In previous work, |Baneriee fc Pudritz|[2006l . l2007l ) have 
successfully implemented the FLASH AMR code to study 
the effects of ideal MHD in the context of star formation. In 
this paper we take the first step and outline the construction, 
testing, and behaviour of a 3D MHD code with adaptive 
mesh refinement (AMR) that employs ambip olar diffusion. 
Our code is an extension of the FLASH code (|Frvxell et alj 
2000), and is thus general enough to treat many if not all 
of the problems mentioned above. Our code also treats the 
full energy equation and is thus able to model effects such as 
drift heating and re-sorting of energy due to the diffusion 
of the field (important as a source of energy dissipation). 
The added advantage of AMR is particularly essential in 
the study of star-forming clumps which we study in some 
detail in our second paper. 

The extension of AD studies to three dimensions to- 
gether with careful treatment of gravitational forces is an 
important step, for several reasons. First, star forming cores 
do not form in isolation in uniform media - they are as- 
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sociated with fila ments presumably created by systems of 
supersonic shocks (I Allen et al.ll2007l ). Second, cores have as- 
socia ted velocity gradien ts that can be interpreted as rota- 
tion |Caselli et alj|2002al lbh and often (but not always) are 
threaded by fields whose mass-to-flux ratio is near unity 
ijHeiles fc Crutcherl [20051. The collapse of rotating objects 
that are threaded by magnetic fields leads to the creation of 
toroidal mag netic fields that play a central role in launch- 
ing outflows dBaneriee fc Pudritj200dlMachida et al.ll2007l ; 



iHennebelle fc Fromand 20081 ). However, toroidal fields and 
rotation have been left out of most simplifying mod- 
els of AD and core evolution, or not fully imple- 
ment ed (Fiedle r fc Mouschovi as 1993; Basu fc Mouschoviasl 



19941; iDesch fc Mouschoviasl l200ll ; ICiolek fc Basul l200rj : 



Tassis fc Mouschoviasll2007l ). Third, on the level of disc evo- 



lution, bars and spiral arms in discs can be very important 
in transporting disc angular momentum - which can be com- 
petitive with extraction of angular momentum by ma gnetic 
tower flows and disc winds l|Baneriee fc Pudrit j|200rj ) . The 
weakened coupling of magnetic fields in the dense interiors of 
discs - where AD will be relatively strong - implies that an- 
gular momentum extraction may only be significant in disc 
surface layers as compared with transport by bars which 
could occur in the bulk of the material. Fourth, fragmenta- 
tion of discs into binaries may be contr olled by the strength 
and orientation of disc m agnetic fields (jBaneriee fc Pudrita 
120061 : iPrice fc Batell2007l) . 

In the subsequent sections we first describe our imple- 
mentation of ambipolar diffusion in the FLASH code under 
the single-fluid approximation, including a brief discussion 
of it's limitations (§2). Next, we outline code tests (§3) that 
we have performed on isothermal and non-isothermal C- 
shocks. Our analysis provides a new important test for am- 
bipolar diffusion codes which accurately account for heating 
and cooling (e.g. using a cooling look-up table instead of 
approximating the equation of state as a step function of 
density), and analytical solutions are developed. We then ex- 
plore the combined effect of gravity, pressure, and magnetic 
fields with AD through a simplified version of the quasi- 
static collapse of a uniform, initially magnetically subcriti- 
cal, gravitating sphere (§4). We compare the results for pure 
hydro, ideal MHD, and our AD MHD and find interesting 
differences in the radial density profiles, which agree well 
with other related published results. We show that our code 
performs very well in difficult test problems, and is there- 
fore an excellent tool with which to explore the full range of 
hydromagnetic processes that are relevant to stellar forma- 
tion, binary formation, and perhaps even the formation of 
massive planetary systems by such processes. 



2 METHODOLOGY 

As a general comment, our approach is to express the equa- 
tions of partially ionized, magnetized, and self-gravitating 
gas is such a way so as to explicitly allow V ■ B 7^ 0. 
We adopted this method because our numerical scheme 
(FLASH) does not obey flux conservation, but instead, op- 
erates by keeping values below truncation-level error. As 
will be discussed in a subsection below ( H2.6p . we found after 
much experimentation that by allowing for non-zero values 
of the divergence of the magnetic field in the equations, and 



then allowing the divergence cleaning procedure within the 
code to remove non-zero values (to below truncation error) 
after each time step, that greater accuracy and numerical 
stability was ensured for our code. 

The potential generality of our code was subjected 
to scrupulous testing. Analytical solutions to problems 
involving ambipolar diffusion are not particularly easy 
given the fact that terms effectively contain second 
order derivatives in the single-fluid approximation (or 
involve two or more interacting fluids in the multi- 
ple fluid picture) . Common tests include the forma- 
tion of stab l e isothermal C-shock s dMac Low et al" I Il995l ; 
iFallei 120031 ; lO'Sullivan fc Downesl 120061 ; iLi et all I2OO6I L 
a nd performing a similar qua si -static collapse to that 
of iFiedler fc Mouschoviasl (119931) dSafier, McKee. fc Stahleil 



ll997l ; lHosking fc Whitworthll2004bF 

The former is an excellent test with analytic solutions, 
although isothermal. We present herein an adaptation of 
the C-shock test which allows for testing of a code with 
ambipolar diffusion energy terms. This differs only slightly, 
though critically^ from the non-isothermal analytical solu- 
tion of (|Wardleifl991al ). The collapse test is a more qual- 
itative test and involves measuring radial profiles and the 
time-scale of collapse when AD is present. We note that we 
have already rigorously tested our MHD module in several 
other papers in the context of magnetized collapse problems 
in both HD and ideal MHD. 



2.1 Theory of ambipolar diffusion 

Consider a gas composed of both ions (subscript i) and neu- 
trals (subscript n). Electrons present in the fluid are con- 
sidered well coupled to the ions in the regime we wish to 
study. The MHD equations for this two-component system, 
in which only the ions couple to the field, are, 



dt 

dpi 
dt 



+ V • (p n u„) = 
+ V • (piUi) = 



(1) 

(2) 



+ ^ " {p n u n u n + P n ) = -p„g - f f (3) 



d(piUt) 



B 2 



+V • p l u i u l + Pi + - BB = -pig 

at \ 2/io po J 

+ f f -—B(V-B), 

po 

(4) 

where p represents density, u is a velocity, P is a pressure, B 
is the magnetic field, g is the gravitational acceleration and 
f f is the frict ional force density from ion-neutral collisions 
|Spitzerlll97"sl ), 



// = lAUPipn K — Ui) — 



1 



PoPao 



where Bad = — and the drift velocity uj is defined 



U d — Ui — U„ 



(6) 



The constant 7ad 



3.28 x 10 13 g" 1 cm 3 s _1 



m 1 -\-m n 
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represents the coupling of the neutrals and ions (the drag 
coefficient). Ions are considered to be typically HCO + or 
Na + which have similar masses (about 29.0 a.m.u.) and 
collision rates with H2 (< av > n i = 1.7 x 10 -9 cm -3 s _1 
|McDaniel fc Masonl [11)731 ')'). The value of 1.4 in /?ad arises 
from the fact that we have about 10% He per H atom in our 
gas dFiedler fc Mouschoviasl Il993l ; iHosking fc Whitworthl 
l2004bT l and is 1.0 in the absence of He (such as in our 
C-shock test runs). Helium is heavier than H2, and thus 
changes the collisional dynamics of the neutral gas. 

Equations {TJ and ([2]) are the continuity equations of 
the neutrals and ions respectively, expressing conservation of 
mass. Equations ([3} and Q are expressions of conservation 
of momentum. Note that while the ions undergo direct MHD 
forces, the neutrals do so only by collisions with the ions 
through the friction term in the HD momentum equation 
©■ 

The induction equation - coupled to the ions - is, 



dB 

dt 



= -V • {mB - Bu,) - (V • B) m 



(7) 



This defines the evolution of the magnetic field. 

Expressing conservation of energy we find the rest of 
our initial equations, 



dE n 

dt 



+ V • [u n {E„ + Pn)] = -U n • f f + p n g ■ U n (8) 



dEj 
dt 



+V 



B \ I 
m ( Ei + Pi + — I + — (m ■ B) B 

2po / po 



+ p i g-u i (v,i • B) (V • B) 

po 



= u * ■ ff 



(9) 



where we see again that the neutrals undergo purely HD 
forces while the forces on the ions are MHD. The energy 
densities are respectively defined as: 



E n — "Zpnlln H 7 Pn 

2 7—1 



Ei = ip^u? + — l —P + 7t— B 2 , 
2 7—1 2po 



(10) 



(11) 



where 7 is the adiabatic index of the gas (we consider equal 
indexes for both gases). 



2.2 Ionization considerations and approximations 

One can simplify these equations by noting that the ioniza- 
tion fraction in molecular clouds and protoplanetary discs is 
very low {pi <C p n ). Subsequent conditions are that the ions 
are perfectly coupled to the field and that the ion-neutral 
collisional timescale (the timescale for an ion to collide with 
any neutral particle) time is m uch less than other physi- 
cal timescales in our system (e.g. iLi et alj|2006l ) so that the 
fluids are sufficiently coupled. This can help effectively elim- 
inate many of the ion equations (see i|2.3f) . Although we will 
always be left with pi in the /3ad term which determines the 
strength of the ambipolar diffusion (lower ion density means 
stronger ambipolar diffusion). 

This approach gives us the advantage of a scheme that 
decreases the number of numerical calculations compared 
to many- fluid schemes, in addition to a relatively simple 



implementation. However, there are limitations to this ap- 
proach. Primarily, there are situations in star forming re- 
gions in which the assumptions break down. These include 
small, dense regions in collapsing cores with densities that 
exceed n n > 10 10 cm -3 , and in certain shocks. When col- 
lisions heat up the flow in a shock so that the neutral flow 
becomes subsonic, the shock may undergo a short, sharp 
transition (sub-shock; iDraine fc McKedll993r ). When this 
happens, the shock is called a J-shock. The ion inertial equa- 
tion is important in the formation of this structure, though 
schemes whi ch neglect ion-inertia have been successful in 
modelling it (|Fallell2003"h . This sub-shock structure will play 
an important role in the chemistry and heating of the shock. 
Strong J-shocks are th ought to be common in star-forming 
regions (|Chernofi1ll9"87h . 

In C-shocks, where the transition of the shock is su- 
personic and everywhere co ntinuous, the flow is suscepti- 
ble to the Wardle instability (|Wardld(l99(jl '). This instability 
is dependent on the continuity equation of the ions, which 
is neglected under the single fluid approximation. The in- 
stability - caused by the buckling of field lines along the 
shock front - leads to the formation of thin sheets (in 3D) 
with density enhancements of about 10 2 times the ambient 
density. These thin sheets - with thickness on the order of 
Iwi = 10~ 2 I/ B hock ~ 2 x 10 13 cm in astrophysical simu- 
lations (|Mac Low fc Sniithlll997T ) - will most likely play a 
role physically in heating the shock and in creating ener- 
getic particles. The sheets are too thin, however, to play a 
role in creating dense cores whi ch co u ld co ntribute to the 
CMF, as has been suggested by iTothl <l 19941). For example 
by looking at the simulations of Mac Low fc Smith] (|l997T ) 
it can be estimated that the Jeans length of the thin sheets 
will be (using their parameters) 



1.7X10 1 



0.01 km s" 1 



\W cm- 3 / 



~l/2 



cm ^> l\ 



(12) 

Thus, the Wardle instability will play an important role in 
mediating the radiation profiles of shocks, but not play as 
important a role for the larger scale dynamics that we wish 
to study. 

The limitation of the single fluid approximation can be 
expressed more ex plicitly for reg ions where the ambipolar 
Reynold's number (Zweibel 2002) 



uL B u 
-Had — — — — — < 1, 



(13) 



where it is a typical velocity and Lb is a typical length- 
scale on which the field varies. Under these circumstances 
the ions and neutrals wil l not be sufficiently coupled to 
allow the approximation (|Padoan et al.l |2000| ) , an impor- 
tant regime for the dissipat ion of astrophysical turbu- 
lence (lOishi fc Mac Low 200tjj). In astrophysical situations, 
it can be shown that u > Ud is generally valid. Fol- 
lowing IZweibel fc Brandenburg! (|l997f ). we assume Lb ~ 
(10 20 cm~ 2 )/n„ based on observed column densities of dif- 
fuse clouds. This allows us to derive a value for Ud/u us- 
ing equation (|14p for the ionization, a typical flow veloc- 
ity of it w 1.5 km s" 1 {M = 5) and a magnetic field 
strength of 10 fiG. The relation satisfies Ud/u < 1 for all 
n n , and for the most part satisfies Ud/u <s£ 1. For example, 
Ud/u = [6 X 10~ 3 ,0.5, 3 x 10~ 2 ,3 x 10~ 4 ] for neutral den- 
sities of n n = [10 2 ,2.1 x 10 3 ,10 6 ,10 1 



|l01 cm 3 , where n n 
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2.1 x 10 3 cm -3 gives the maximum value of (itdAOmax = 0.5. 
Thus the approximation is a valid one on the scales we are 
interested in. 

In this paper, we shall use a simplified approach in treat- 
ing the ionization of molecular clouds which is pertinent to 
lower density regimes, and which has been used widely in the 
literature of the subject. The ion density can be expressed as 
a simple function of neutral d ensity in this regime. We fol- 
low the lead of other authors (jFiedler fc Mouschoviaslll993l : 
iHosking fc Whitworthl l2004al ) and use the very simple ex- 
pression (see the quasi-static collapse problem): 



= K 



( n n \ fc „, t n n y 
V10 5 cm-3/ ^ Vl0 3 cm-3y 



(14) 



where n is a number density, K — 3 x 10 cm" , k — = 
and K' — 4.64 x 10 -4 cm -3 . This approximation is meant 
for density situations in which dust grains do not play a 
substantial role in the ionization balance. It allows one to 
eliminate entirely the need to track the ion density in the 
single fluid approximation. The second term dies off quickly 
in the higher density regime, at which point we're left with 
the common nt oc nl/ 2 relation. One can, by these means, 
reduce the 2-fluid (ions plus neutrals) to that of an effec- 
tively single fluid with low conductivity. We note that this 
approach is still very important and useful since it allows 
one to accurately track the physics of core formation and 
collapse up to densities that far exceed the density of cores 
that are nearly in equilibrium (at 10 4 — 10 6 cm" 3 ), which is 
our intention. 

This well-known expression arises by app roximating 
the r e sults of ionizat ion equilibrium calculations (|Elmegreenl 
Il979l ; iNakanol [l979) where the sole form of ionization is 
through cosmic rays (where the cited authors assume a 
cosmic ray ionization rate of fo = 6.9 x 10~ 17 s" 1 ). It 
is approximately correct for densities n < 10 10 cm -3 , af- 
ter w hich charged grains become the principal charge car- 
riers (Tas sis fc Mouschoviasl l2007f) . Also, more active and 
complicated chemistry will be required to model C-shocks 
(|Pineau des Forets et alj|l986h . However, for our purposes, 
it is sufficient that it matches th e equation originally used 
bv lFiedler fc Mouschoviasl i|l993h . Furthermore, in order to 
study the collapse and disc formation to densities greater 
than 10 10 cm -3 a multi-fluid approach including grains will 
be required. 

For the C-shock test we use a constant ion density. This 
leads to slightly different analytic solutions than used in 
the past. We demonstrate how to adapt given a certain ion 
density prescription, or non at all (ion mass conservation). 



2.3 Carrying out the single fluid approximation 

Our goal in this Section is to make equations ©, ©i 
and © look like MHD equations for the neutrals plus some 
other terms. This allows us to simplify the computations 
significantly by reducing the number of equations that must 
be solved at each time stejQ . We have already noted that 



More general and robu st mult i -nuid techniques have been 
developed llO'Sullivan fc Downes] |2007|. |2006|; iLi et al.l l2006t 
lOishi fc Mac Lowl2006l ; lFallel2003l) , although isothermal and on a 
static grid. These techniques may be more important in turbulent 
conditions where pertinent timescales are highly variable. 



molecular clouds are dense and poorly ionized which implies 
that pi <C p n . Furthermore, we require that the region we 
study have an ion-neutral collisional time (tm = l/(7ADPn)) 
that is much shorter th an other physical timescales in the 
problem (|Li et al.ll2006l ) . This is also typically true in molec- 
ular clouds, which allows us to ignore the ion inertial forces 
with respect to their magnetic and frictional forces. As is 
well known, one can use this to reduce the ion equations of 
motion into one equation for the drift velocity, saving half 
the calculations. 

More explicitly, we find that the gravitational and pres- 
sure forces for the ions are negligible with respect to the 
magnetic force f m — J X B (where J = ^j-V X B). We 
can then ignore the ion inertia through our secondary ap- 
proximation. Equation ((4} leaves us the equation of force 
balance - between the full Lorentz force exerted directly 
upon the ions, and the friction force that arises from ion- 
neutral collisions: 

ff = -fm = ~J x B, (15) 
from which we can solve for the velocity of the ions: 

= /3ad(V X B) x B (16) 



ut =u n + /3 A d(V X B) x B. 



(17) 



Substituting (1171) into © we can solve to get [simplify- 
ing V ■ (uB - Bu) = -V X (u X B)], 

— =V x (u n x B) - (V • B) u n 

+ V x {Pad [(V X B) x B] x B} ( 18 ) 
- (V • B) [/? AD (V X B) x B] , 

where the neutrals look to be perfectly coupled to the mag- 
netic field save for the two last terms which describe the 
ambipolar diffusion of the field. 

Similarly, in the energy equation of the neutrals © we 
substitute into the frictional force (11511. We recover, 



dE n 
dt 



+ V 

+ V 

1 

Mo 



B 2 \ 1 
E n + P n + — ] + — (u n ■ B) B 



PadB' (J x B 



2mo / Mo 



(u n • B) (V • B) + po/3ah || J x B\\ 



Pad [B-(Jx B)] (V • B) 



(19) 



where the last term on the LHS and the last two terms of the 
RHS are ambipolar diffusion terms. Heating of the fluid oc- 
curs through the dissipation of the field (/^o/3ad || J X i?|| 2 
term). Diffusion of energy is similar to that of Ohmic dif- 
fusion (V • [/3adB 2 {J X B)] term), given the ambipolar 
diffusivity below. The rest (ignoring V • B terms) looks like 
the MHD energy equation for the neutral species; the total 
energy of the fluid is changed by the magnetic field as it 
undergoes motion in the gravitational field. 

Now that the single fluid equations have been written 
out, we may dispense with the subscript n (i.e. representing 
the gas as a whole). The following MHD equations are now 
in a form that is optimized for our computational approach 
(i.e. similar to IPowell et al.l |l999)), wherein fluxes appear 
within a divergence on the LHS and sources are placed on 
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the RHS: field to form which lead to instabilities generated by the am- 

Qp bipolar diffusion terms. This creates unphysical states such 

— + V • (pu) — (20) as ne g a ti V e densities. 



d(pu) 
dt 



B 1 

puu + P H BB ] = -pg 

2/i po 



— B(V B) 

Po 



2.5 The FLASH2.5 code 



dB 
~dt 



+ V • (uB - Bu) 

+ V • {poPad [(J X B) B — B (J x B)]} 
= - (V • B) u - (V • B) [poPad(J x B)] 



(21) 



(22) 



8E 
dt 



+ V 



u[E + P + 



2p 



Po 



(u • B) B 



+ V • [PadB 2 (JxB)]=pg-u-^-{u-B)(V- B) 

+ Po/3ad II J X B\\ 



pit 

'- ■ /3ad [S • ( J x B)] (V • B) 



(23) 



2.4 Diffusion and ambipolar drift: typical 
timesteps 

A glance at the form of the induction equation above (|22p 
clearly shows that ambipolar diffusion (i on-neutral drift) is 
not just the simple diffusion of the field llOishi fc Mac Low! 
120061 ; iMac Low et~ail (l995l ; iBrandenburg fc Zweibellll994 ). 
Ohmic diffusion, for instance, only adds a term of the form 
V x (r;V X B) to the LHS of the induction equation (under 
a similar approximation), where r/ is the Ohmic di ffusivity. 
The c orresponding ambipolar diffusivity is given as (jZweibell 
[20(1), 



IJAD = PadB , 



(24) 



where jjad has units of cm s . This can be used 
to re-evaluate the ambipolar diffusion induction term as 
l|Brandenburg fc Zweibeilll994r i. 



V x [n AD V X B po/3av>(J • B)B] 



(25) 



on the LHS. Clearly the second term demonstrates the devi- 
ation of ambipolar diffusion (ion-neutral drift) from a purely 
diffusive process. 

In fact, it is easy to imagine situations in which J ■ B 7^ 
0, such as disc accretion in a collapsing core. It would then 
not be surprising that the extra non-diffusive term becomes 
important as field lines are strongly oriented against the ac- 
cretion flow, which caries a charged current. Nonetheless, 
we expect the diffusive component to dominate. We thus 
use a diffusive scaling to estimate the typical time-scale 
that is needed in ord er to satisfy the Courant condition 
l|Mac Low et al.lll995h . 



TAD ~ To 



(Aa 



I)AD 



(26) 



where the factor of To is a fudge factor (we use 

Tp = 0.00833 « 1/120 in our simulations as used by 

IMac Low et~ail (|l995ln and Ax is a typical length scale (see 
^A3|) . Longer time-scales or slower speeds (if the timestep is 
not used, for instance) tend to allow sharp gradients in the 



The FLASH code l|Frvxell et al.ll2000h employs a wide range 
of HD and MHD physics. It has been well tested and has 
become quite popular in the astrophysical community. Its 
principal advantage is that it can perform AMR, as this 
saves considerable computational time while maintaining a 
high level of refinement. A secondary, though important ad- 
vantage for astrophysicists, is that it does a very good job 
in capturing shocks. 

The refinement criteria fo r FLASH have been adjusted 
somew hat in previous work by Baneriee . Pudritz. fc Holmes! 
(2004). It was shown bv lTruelove et alj (J1997) that one needs 
to refine the Jeans length, 



A, 



Gp 



1/2 



(27) 



by at least four cells in an AMR code, where c s is the isother- 
mal sound speed, G is Newton's constant and p is the lo- 
cal density. If this is not satisfied artificial fragmentation 
induce d by the numerical grid will result. iBaneriee et al.l 
l|2004l ) setup the refinement criteria such that a variable frac- 
tion of cells will resolve the Jeans length. In the quasi-static 
collapse we use 8 cells per Jeans length, guaranteeing artifi- 
cial fragmentation will be avoided. 

Our explicit coding of the above single-fluid approxi- 
mation in FLASH is described in detail in Appendix |Qy. 
This includes the splitting of the equations, the timestep 
condition and the dimensionality of the equations. 



2.6 V • B considerations 

The MHD module in FLASH employs a somewhat contro- 
versial sc heme that do e s not explicitly constrain V • B = 0, 
based on lPowell et""aH (119991 ^1. The V ■ B terms are left in 
all of our equations to emphasize this fact: V • B is not 
explicitly zero in FLASH, but c onstrained below tr uncation- 
level error after each timestep l|Powell et al.|[l999r ). Consid- 
ering these terms during a timestep is thus important in 
guaranteeing the stability of any computational scheme ap- 
plied to FLASH'S MHD module. We find it helps decrease 
V • B by a couple orders of magnitude in longer non-steady 
computational runs, however it may not be important in 
more steady flows. Note that we are also using a method 
employed in FLASH which eliminates V • B as it's created 
through a simple diffusive method (this is step is performed 
after each timestep). 

The evolution of V • B in our isothermal C-shock test 
(33| evolves very much the same with or without these extra 
terms (there is only a minor advantage of having such terms 
in the code). As the shock reaches a steady state, the aver- 
age value remains constant. We note that this former case 
involves small integration times of a few hours and a steady 
final state. Our simulations of the pre-stellar collapse of a 

2 The code has been well tested and, despite the controversy, is 
widely used in numerical AMR— based MHD simulations. 
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Bonnor-Ebert sphere performed in a companion paper is 
continually evolving over integration times of a few days. It 
shows orders of magnitude improvement when these terms 
are added and more than 10 orders of magnitude improve- 
ment to similar runs performe d by the SPH simulations of 
iHosking fc Whitworthl (|2004bh . In a general sense it may be 
advisable to only use these terms if V ■ B becomes an issue 
(e.g. non-steady states, long integration times). Otherwise, 
skipping the lines of code pertaining to these terms could 
save important computing time. 



3 TEST CASES - OBLIQUE C-SHOCKS 

To test the code we have performed a detailed check with 
an isothermal C-shock in addition to a second C-shock test 
which we believe is the first quantitative test in practice 
of an ambipolar diffusion code with a full non-isothermal 
energy equation. Our hope is that this latter test is applied 
to future codes developing ambipolar diffusion with proper 
energy considerations. 

C-shocks are continuous transitions of unconserved HD 
variables mediated by elastic collisions be tween charged 
and neutral particles ( Draine fc Mc Kcc 1993). Without am- 
bipolar diffusion (e.g. ideal coupling) such flow transitions 
are highly discontinuous. C-shocks have quickly become 
a common test case for ambipolar diffusion codes due to 
recent interest that has sprung up in a variety of sub- 
jects , particularly in turbulenc e models |Mac Low et al.l 
19951 ; ISmith fc Mac Low! Il997l ; iMac Low fc Smith! Il997t 



Falldl 20031 ; lO'Sullivan fc Downedl2006l ; iLi et alj|200^ . Tests 
for ambipolar diffusion, with a given analytic solution, are 
very hard to generate due to the complex nature of the equa- 
tions. The C-shock has provided a very convenient test case 
for the papers mentioned, all of which are isothermal. As as- 
trophysical problems generally involve shocks, we will even- 
tually require full energy considerations with heating and 
cooling. A test for such simulations which include the ef- 
fects of drift heating and energy diffusion does not exist to 
our knowledge. However, further development of the idea 
of C-shocks can be expanded to include these effects. We 
present the first such test of a non-isothermal ambipolar 
diffusion energy code that we are aware of in the literature, 
in the hopes it will be used as a standard test in the future. 
Alongside this we perform an isothermal C-shock test to 
independently test the induction equation. 

3.1 Basic setup 

In Fig. [T]the basic idea of an oblique C-shock is illustrated 
in our initial setup, such that the field B is initially oblique 
to the normal of the shock front by an angle 9 S ■ The disconti- 
nuity will be transformed into a more continuous transition 
as the neutrals sift through the ions and feel the magnetic 
forces through collisions (modeled as friction) . A final steady 
and continuous transition is obtained. The dimensionality of 
the problem is two, so we run the simulation in permutated 
orientations of x, y and z to fully test our code. Nonetheless, 
each run involves a 3D tube. The boundary conditions (for 
the standard permutation presented here) are inflow at the 
left x boundary and outflow at the right x boundary, while 
periodic conditions are used for every other boundary. For 



ID plots we take data from a central row of cells in our box 
(though no particular row appears different). 

The C-shock ana lytical solutions a r e based on work 
from the 80's and 90's (lDraindll980l , Il986l ; IWardle fc Draine! 
ll987l ; IWardlell990l . ll991bl lal). where the first paper deals with 
the oblique shocks we present here. The equations presented 
therein offer a working analytical solution to a steady state, 
isothermal, two-fluid C-shock, mediated by ambipolar dif- 
fusion. 

In respect to the equations we present below, we have 
found some inconsistencies with how the energy equation 
was reduced in the past to give a differential equation for the 
pressure. Diffusive contri butions to th e MHD energy equa- 
tion were neglected in llDrainel Il986l) . altho ugh frictional 
heating can be added as in IWardld dl991al) ( in addition 
to cooling). Note however, IWardle fc Draind i|l987l ) argue 
that the electromagnetic field cannot transport energy as 
E X H = 0. With this in mind, they find the magnetic 
terms in the energy equation (save for drift heating) will 
cancel out when the single fluid approximation is carried 
through and the shock conditions are applied (in the E = 
frame). We show below that this is not the case when dif- 
fusive terms are considered and that magnetic terms indeed 
play important roles in the energy equation. The agreement 
of our test results with our theoretical predictions and the 
significant differenc e between our p redicted pressure distri- 
butions with that of IWardle] (|l991af ) certifies our viewpoint. 

Also note that since we use a single fluid code, we must 
make a choice on how to treat the ion density evolution. A 
simple choice would be to keep it constant (p; = p; ). To 
account for this involves adjusting a term in the analytical 
solution, and we outline how this is done below (as well as 
how to change the equations to ion mass-conserving form 
or something else). 

We take the initial conditions of IMac Low et all (|l995l ) 
as their shocks show long transitions of (10 — 20) Z/ s hock, in- 
dicating strong ambipolar diffusion. Running these through 
the analytic solutions we find the values in the post-shock 
region. Our initial setup consists of a left state with these 
pre-shock values and a right state with the post-shock val- 
ues (Fig.[T]and Table[T}. We evolve the shock tube and allow 
it to settle to its final configuration, which we compare to 
analytical predictions. 

The analytic solutions we derived, found by solving the 
steady multi-fluid equations for the final shock, consist of a 
coupled set of two first order ordinary differential equations 
in parameters p and b. The e quation s are, w here we retain 
the parameter conventions of lWardld l|l991al ): 



- "/r n p 



(7- 1> 

7ADPi 



dp 

dz 



s n + sin ( 



r + 



(G» " An) 



db = 7ADPi A 2 
dz v s 



(28) 



(29) 



where, 
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Figure 1. Our initial setup of the oblique C-shock. The left, pre-shock state is 4 -L s hock in length while the right, post— shock state is 
8 Lsjjock i n length. The y and z lengths are both (1—2) I/ s hock i n length, for high and low resolution runs respectively. After a while the 
shock will smooth out the discontinuity and settle to form a static C-shock. 



Table 1. The initial conditions for the oblique C-shock test. The value of 7ad is 1-0, and the ion density is kept constant at pi = 
10~ 5 cm -3 . The temperature is adjusted to 10 K by setting the molecular mass appropriately, given that the sound speed is c s = 
0.1 cm s . The initial field deflection is 9 a = it/ 4 and B x = Bo cos (9 3 ) = v4?r cos (9 a ) is a constant. Non— isothermal post-shock states 
are only slightly different than isothermal ones (save pressure). 



state 


Pn [g cm 3 ] 


P n [dyne cm 2 ] 


v x [cm s 1 ] 


v y [cm s 1 ] 


By [G] 


left 


1.000 


0.0100 


v s = 5.000 


0.000 


2.507 


right iso 


8.045 


0.0804 


0.621 


0.840 


23.553 


rlght non _ iso 


7.976 


0.5 


0.627 


0.830 


23.313 



l-(P-Po)-P#) 



r = I 1- — 



b 2 + cos 2 



br„ (s n + sin S ) + cos S 



bo 2 
cos ( 



A 2 



(30) 



(31) 



(32) 



(33) 



and 7 is the adiabatic index of the neutral gas, 7ad is the 
collisional constant between ions and neutrals described in 
Sj2T] p no and pi are the initial density states of the neu- 
trals and ions respectively, v 3 is the initial speed of the gas 
towards the shock, A — v 3 /va is the Alfven number of the 
pre-shock state and va = Bo /^/47rp„ is the Alfven velocity. 
The velocity compression of the neutral and ionized gases are 
described by the parameters r,: = Vs/vi x and r n = v s /vn x . 
The dimensionless pressure is p = P n j (pn fs), and the di- 
mensionless field (in the y-direction) is b = B v /Bq, where 



Bo is the total initial field magnitude. The deflection of the 
field is quantified by s n — (v ny B x ^) / (v s Bo), as the deflected 
field will accelerate a velocity in the y direction as it moves 
through the gas. The x-component of the field stays con- 
stant due to the symmetry in the problem, and in our runs 
the ionization is kept constant (though ion velocity shocks 
strongly). The parameter r defines what type of chemistry 
is used or if ion mass is simply conserved. The value given in 
Equation (|31[) is for a constant i on density. For ion mass con- 
servation (as in I Wardkl (|l991al )) the value is r = (n — r n )- 
For a general chemistry, r = f(x) (1 — r n /ri) (such that one 
can derive pi — f(x)pi , where f(x) is some function in the 
x direction which may include variables such as r n (x)). 
Drift heating is governed by the parameter: 



G n = JADPiPn \\ 



v, 



= 7AD PnoPioV 2 ,-^- (b 



PoPao II J X B\ 



2 + cos 9 3 2 ) , 



(34) 



and A n is the cooling rate. The analytical solutions prove 
to be rather difficult to solve without cooling (via a simple 
Runge- Kutta 4 techniq ue) . By including a cooling rate as 
done in IWardld (|l991al ). solutions can be obtained that are 
significantly different than un-heated shocks yet realistically 
integratable. We take a simplified cooling rate that approx- 
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imates those of iLepp fc Shull (|l983l ), relevant to molecular 
clouds (though meant purely to simplify the computation of 
this analytic model): 







[erg s 



(P > CcoolPo) 

(otherwise), 



(35) 

where po is the pre-shock dimensionless pressure, e coo i de- 
fines at what pressure the cooling turns on and A no is a 
fudge factor. We found A„ = 5 x 10 and e coo i = 50 pro- 
vide good analytical graphs at a reasonable computational 
cost, while encapsulating the shock in a similar lengthscale 
to the isothermal version of the C-shock. It was important 
in our case to keep the shock length small as we ran the 
shock tests in 3D cartesian coordinates. 

The lengthscale of the shock is given as L s hock = 
\/2fMifiow, where ifl ow = 1/(7ad/9j ) is the time-scale for 
the shock. We use these parameters in our simulation setup, 
having a box 12 x (1-2) x (1-2) Z/ a hock in dimensions (depend- 
ing on resolution), and printing plot files every 0.1 tfl ow (for 
a total of 500 files). The evolution consists of a large outward 
propagating density enhancement which is ejected from the 
box as the flow settles to form a C-shock (we consider the 
C-shock 'steady' at about 20 ifl ow )- Note that while this 
density wave is in the box, timesteps are smaller by about 
an order of magnitude. Because of this, a longer box means 
longer computation times. 

Note that our solution for the pressure derivative in 
equation (I28|) is significan tly diffe r ent th an that derived in 
iDrainel i|l986n and used in lWardld (|l991al ). The difference is 
the additional term on the right hand side of Equation (128[) : 



s n + sin ( 



(36) 



We compare our solution to that of I War die] (|l991aT ) in 
Fig.[2]for an adiabatic index of 7 = 1.4 and the same cooling 
rate defined above and using the values presented in Table 
[T] The significant difference between the results is quickly 
apparent. See typical isothermal values of pressure in Table 
{TJ or errors in Fig.[4]to understand the scope of the numbers 
in this figure. Other quantities have only minor differences. 

3.2 Numerical results 

In the following subsection we present comparisons of 
isothermal and non-isothermal C-shock simulations rela- 
tive to analytical sol utions. The error is calculated as in 
iMac Low etaD (|l995h : 



% error = 100 



analytic ^numerical 



(37) 



max(g ana i ytic ) 

We employ rather low resolution in these simulations; our 
'low' resolution has cell sizes of (1/2) L s hock while our 'high' 
resolution run resolves (1/4) L s hock (though in 3D). For all 
our models we plot the percent error for the high and low 
resolution cases. It can be seen from the figures that the high 
resolution error is strongly reduced in relation to the low res- 
olution errors. This is enough to show the high accuracy of 
our m odels (giving similar errors to that of IMac Low et all 
|l995h with half the resolution) and demonstrate a clear 
convergence, thus confirming our implementation of the in- 
duction equation. Our results are in Fig. @- 



Table 2. The parameters used in the quasi-static collapse test 
run. Note that the ionization is the same described in i|2.2l where 
/i mo l = 2.33. The sound speed then corresponds to a temperature 
of 10 K. The box dimensions extend from —3.24 pc to +3.24 pc. 



-Rtest [pc] n in [cm 3 ] riout [cm 3 ] b tC st [/^G] c s [km s l ] 



0.75 



300.0 



30.0 



30.0 



1.889 



Finally we present the results of the non-isothermal 
run, testing our energy terms. For these simulations we use 
an adiabatic index of 7 — 1.4. These simulations also show 
high accuracy and a clear convergence. Our results are seen 
in Fig. @ and @. The sharp edge in the post-shock region 
of Fig. is a result of adjusting the parameters in our toy 
cooling model (|35|l . This was done to properly encapsulate a 
stable non-isothermal C-shock on a grid size similar to the 
isothermal C-shock. 

It is clear from the above tests that our code performs 
admirably in all facets. 



4 TEST CASE - QUASI-STATIC COLLAPSE 

This test involves matching qualitatively the characteristics 
(such as time-scales and curve features) of the quasi-static 
collapse of a thermally critical, but otherwise magnetically 
subcritical, non-rotating core. In or der to study the effects 
of AD i n an ax isymmetric geometry. iFiedler fc Mouschoviasl 
<| 19931 ) (|FM93l hereafter) set up an initial state with a peri- 
odic boundary conditions (along all surfaces) meant to rep- 
resent an infinite chain of fragments. They used a static, 
specially designed adaptive grid to follows the evolution. Its 
history as use for a test code for amb ipolar diffusion is spo- 
radic, appearing in lSafier et al.l ljl997h (an analytical model) 
and mo re recently in SPH models bv lHosking fc Whitwortbl 
i|2004bh . 

Our purpose is to investigate a test problem whose col- 
lapse will have similarities with that of the well-studied 
Bonner-Ebert models. There is consi derable evidence that 
these latter states actually exist ()Alves. Lada. fc Ladal 
2001). Nevertheless, in order to make comparisons with the 
IFM93I calculations, we use the physical parameters from 
their model 1 within our own spherical initial states. We 
present these physical paramete rs in Ta ble [2] Although the 
model we use is different from FM93, it produces a thin 
quasi-static disc of similar characteristics. In this sense it 
is studying the same problem, the time-scale between thin 
disc formation and dynamic collapse. The difference in mod- 
els is due to th e uniqu e way in which the problem was ini- 
tially setup in IFM93I . using periodic boundary conditions 
on a small cylinder - whereas we employ an adaptive mesh 
refinement scheme in a 3D Cartesian grid. In the end our 
model studies the problem of an isolated cloud, while the 
cited authors studied closely packed clouds. 

Our model consists of a sphere of uniform density ni n = 
300 cm" 3 , radius R tCBt = 0.75 pc w 1.6 x 10 5 AU , and total 
mass 30.3 Mq. The critical mass for th e inner sphere if it 
were non-magnetized is M cr = 6.6Mq (|McKee fc Ostrikerl 
2007), so the sphere is thermally unstable. However, the 
initial mass-to-flux ratio is V = 0.298 < 1, meaning the 
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Figure 2. Pressure profiles comparing our analytic solution to that of iDrainel [|l986h and IWardlel dl991al) . A significant difference is 
apparent. 



cloud is supported by its magnetic field. The sphere is in 
thermal equilibrium with an external medium of lower den- 
sity (pambicnt = 0-lpcorc). The s etup is shown schematically 
in Fig. © The model of IFM93I used a uniform cylinder of 
density p C orc, radius J?tost and height above the midplane 
Z = -Rtost- Both models have periodic boundary conditions, 
thus the latter will have less of an ambient source of material 
in which to accrete due to it modeling a closely knit cluster 
of collapsing thin discs (as a consequence of a small box size 
and periodic boundary conditions). Also, our initial setup 
will collapse much quicker into the qu asi-stat ic state as it is 
initially highly unstable. The setup of IFM93I is only slightly 
thermally critical to fragmentation. We run pure HD and 
ideal MHD models alongside the non-ideal MHD model to 
shed light on what is happening physically. 

The sphere quickly collapses along the field lines into 
a quasi-static disc (which matches the description of the 
thin disc model described by IFM93I and many others) . The 
evolution of the central density is shown in Fig. [7] The initial 
rise in density i s due to the rebounding after the initial infall, 
as described bv lFM93l . Gas piles into the disc while magnetic 
pressure is pushing it out. At this time (about the free fall 
time, tff = 1.96 Myr), without the magnetic field the cloud 
would have entered a freefall collapse (illustrated by the pure 
HD curve in Fig. [7^) ■ However, the disc reaches the quasi- 
stable equilibrium - where gravitational forces are balanced 
by magnetic forces. The infalling gas on larger scales takes a 
while to respond and continues to accrete more matter on to 
the disc. Magnetic pressure forces then drive out the excess 
matter and the disc settles into a quasi-equilibrium. Note 
that the cited authors reach this characteristic 'rebound' 
peak at 6 Myr, while we reach it at 2 Myr. 

The ideally magnetized case will evolve only marginally 
after this, accreting material from the ambient medium. 
This is because the mass to flux ratio of this state is still 
subcritical (The maximal V for the ideal model is about 
r = 0.41 < 1 at the end of the simulation, while it is about 
4.90 by the end of the non-ideal run). The non-ideal disc 
also continues to accrete mass internally as well as from the 
ambient medium - its central density evolves at a higher rate 
than that of the ideal case. This is illustrated by the evolu- 
tion of the density as well as the radial velocity as a func- 
tion of distance, and for different times in the simulations 



(Fig. [SjQ- Of note is the considerable difference in density 
ranges for each model. Inside of 10 4 AU the ideal model has 
very small inflow motions. While, in the later stages of its 
evolution, the non-ideal model appears more like the purely 
HD model, both having considerable inflow motions inside 
of 10 3 AU. This indicates that indeed significant inflow is 
occurring inside a region of the disc that is supported in the 
ideal limit. It is interesting to note how the infall motions 
resemble the outside-i n collapse motions that ch aracterize a 
Bonnor-Ebert sphere (|Baneriee fc P udritz 2 007h . 

The field lines are found to remain straight (until full 
freefall collapse is reached) despite advanced infall in the 
AD model. This corresponds to neutral gas moving, under 
the influence of gravity, through the ions attached to the 
magnetic field. 

Eventually, freefall occurs in the non-ideal model. This 
is because ongoing AD in the quasi-equilibrium state grad- 
ually increases the mass to flux ratio until a cri tical st ate for 
collapse is achieved, and then surpassed. The IFM9&1 model 
has a freefall collapse that begins at about 16 Myr while ours 
begins at about 12 Myr. The difference is readily understood 
by observing that each freefall state happens about 10 Myr 
after the rebounding 'bump' occurs and the quasi-static disc 
is formed. In this sense our models agree on the time-scale 
of ambipolar diffusion (although, by having a larger reser- 
voir of material to accrete from, we slightly shortens our 
time-scale from 10 Myr). 

The radial profiles of Fig. [HJ; and Fig. [SJI also agree 
well with the previous work. FM93 find dlogn/dlogr ~ —2 
while we find a va lue of -1.89. In areas more supported by 
the magnetic field. IFM93I found a profile that goes more like 
-1.6. In our ideal run, we found a profile of -1.51 (the cited 
authors do not list detailed properties of the ideal collapse). 
It is interesting to note that our HD freefall collapse had 
a profile of -2.3, indicating a restricted build-up of density 
in the non-ideal model, clearly due to the magneti c field. 
The radial velocity profiles are also similar to IFM931 . show- 
ing similar ranges in v r and a similar profile. Supersonic 



3 The profiles are radially binned spherical averages, such that 
for a quantity v, v(r) = f v(x)dQ, where dQ = dcos8d<j). 
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Figure 3. Profiles comparing simulations of isothermal C— shocks with analytical solutions at different resolutions for (a) density, (b) 
y-component of magnetic field, (c) x-component of velocity and (d) y-component of velocity. Accuracy is high and convergence is evident. 
Pressure distributions are identical (though scaled) to density. 
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Figure 4. Pressure profiles comparing simulations of non-isothermal C-shocks with analytical solutions at different resolutions. Accuracy 
is high and convergence is evident. Temperature distributions are similar (though scaled) to pressure. 
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T = 10 K 



Figure 6. Setup of the the quasi— static collapse test problem. Our setup of the quasi— static collapse is different than the one presented 
by IFM93 . The key parameters are shown. The whole box is kept at 10 K throughout the collapse. 
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Time [Myr] n c [cm -3 ] 

Figure 7. Results from the first test run showing the evolution density and magnetic field initiated from our quasi-static collapse setup. 
We plot (a) central density evolution and (b) the evolution of the central magnetic field. The ambipolar diffusion run is shown in the 
solid line and the ideal run is shown in the long dashed line. A non-magnetized 'hydro' run is shown as a short dashed line is a). 



velocities are achieved in all cases. We do not detect a shock 
forming in any of our models. 

The evolution of the central magnetic field with density 
is shown in Fig. [7p. Our curve has the right qualitative fea- 
tures: a roughly constant field followed by a positive power 
law. The p ower la w begins at roughly (2-4) x 10 4 in agree- 
ment with IFM9&I At greater densities, there is a distinct 
power-law relation between these quantities - we find 

B c = 0.51n*; « ^ 0.37 (38) 

which compares well, though not exactly, with their result of 
k — 0.47. We note however that we chose to resolve our runs 
by every 0.1 Myr and were left with only two data points to 
calculate k. 

In Fig.|5]we find that the ideal disc tends to be thicker, 
in the sense that density distributions are more spread out in 
the disc. The non-ideal disc has a significantly higher central 
density with a corresponding higher degree of refinement 
(the ideal case has 3 levels of refinement while the non-ideal 
case has 8 by the time we end the simulation at 11.9 and 
11.4 Myr respectively). This is because the non-ideal disc is 
more collapsed due to a slight infall motion present of about 
0.1 km s _1 . Since free-fall collapse has begun in the non- 
ideal case by our final plot file we can observe that field lines 
are being pinched in by radially accreting material (shown 
in Fig. UOp . As we only output plot files every 0.1 Myr, we 
did not follow the collapse further (also cooling restrictions 
as well as shock and drift heating would become increasingly 
important in further stages of the clouds evolution) . 

Also of note is the preservation of axisymmetry in the 
simulation, but also the non-linear outline of the thin-disc. 
The latter is due in part to the collapsing sphere which falls 
more quickly in the centre regions than its sides. This leaves 
peaks and ripples in the disc's surface as the outer parts of 
the sphere eventually join the disc. The disc is also much 
larger in radius than the initial sphere, aiding our argument 
that box-size is important in this type of calculation. 

We ha ve corr ectly found the time-scale of collapse de- 
scribed by IFM93I . matching all qualitative features of the 
quasi-static collapse. We conclude that we have passed this 
qualitative test. More importantly, we hope that we have 



set up a version of the collapse that is easier to reproduce. 
However, optimization is required in the parameters in order 
for this test case to be streamlined for quicker testing. 



5 CONCLUSIONS AND DISCUSSION 

Our results present a properly implemented and highly ac- 
curate application of ambipolar diffusion into the FLASH 
AMR MHD code. Our model is that of a single-fluid (which 
provides a computational advantage over multi-fluid meth- 
ods with a necessary chemical model at the expense of more 
stringent physical restrictions) . We have the benefit of a full 
non-isothermal energy equation. 

We have presented the first test case of non-isothermal 
aspects of ambipolar diffusion through the development of 
the corresponding C-shock analytical solution, with appli- 
cations to a general chemistry (or no chemistry at all). We 
hope this becomes an important test for future codes. 

The quasi-static collapse simulation was modified in 
a more computationally accessible manner. The added ef- 
fect of a larger mass reservoir speeds up the creation of the 
quasi-static thin disc, and slightly shortens the magneti- 
cally mediated time-scale for collapse (by about 0.5 Myr). 
It may also explain the slight differences we found in power 
law trends (power laws tended to be lower by about 0.1). 
Specifically, we found that the uniform spheres quickly adopt 
flat-topped, B-E like density profiles with n(r) oc r~ a where 
a = (1.51, 1.89, 2.30) for our ideal MHD, AD MHD, and pure 
hydro cases respectively. The ideal MHD sphere remained 
supported by the magnetic field (it was initially subcritical, 
although with a more concentrated density profile). Both 
the AD and hydro cases went into collapse. We also found 
the very interesting scaling of the core magnetic field with 
the maximum density B c oc n° 37 , which is in good agree- 
ment with other published results. We found similar qual- 
itative features not only to the comparison study, but also 
to st udies involving the collap se of Bonner-Ebert Spheres 
(e.g. iBaneriee fc Pudritj 120071 ). Although only meant as a 
qualitative test, our quasi-static collapse (alongside previ- 
ously established work) may prove a useful comparison for 
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Figure 8. Radial profiles in the quasi-static collapse test at different times for ideal, non-ideal and hydro models. Profiles are of density 
and radial velocity in the a-b) ideal, c-d) non-ideal and e-f ) hydro runs respectively. Curves are chosen such that they are evenly spaced. 
Values of density and radial velocity are spherically averaged for a given value of r. 



observations and other types of collapse simulations looking 
to test persistence of such a process in our current under- 
standing of star formation. 



In our next paper, we use our new code to investigate 
the collapse of a rotating, magnetized, Bonnor-Ebert sphere 
in full 3D and with AD included. 
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APPENDIX A: IMPLEMENTATION OF NEW 
PHYSICS INTO FLASH 

Here we describe how we have coded our model of ambipolar 
diffusion into FLASH. 



Al Implementation of ambipolar diffusion 

To put the equations in code form we will have to first 
separate sources and fluxes, and then add these to the 
Flux and SOURCE arrays during the sweeping process 
of mhd_SWEEP.f90 in FLASH2.5. This module performs 
the sweeps, calls the interpolation of the data for ideal 
terms and calls fluxes and sources for a given sweep direc- 
tion (including non-ideal data). After this, it evolves the 
MHD equations. In coding the ambipolar diffusion terms 
we follow closely how the resistive fluxes are applied (see 
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Figure 10. Zoomed— in density pseudo-colour from the end— state of the non-ideal quasi-static collapse showing the bending of the 
magnetic field lines. Free-fall collapse has begun and accreting gas in the disc pinches the field inwards. Typical velocity vectors are 
shown in the bottom left of the image. 



mhd_add_resistive_fluxes.f90 in FLASH2.5). A central 
difference method is used to discretize derivatives (V X B 
terms) at cell centres. This should be an effective method 
as the ambipolar diffusion terms, like the Ohmic dissipa- 
tion terms, are parabolic in nature. They won't require the 
highly involved interpolation the rest of the MHD equations 
undergo to be accurate (the ideal MHD equations are hy- 
perbolic). 

To illustrate the central differencing technique, we eval- 
uate an imaginary flux term of the form 



Thus our final flux term will have the form: 



F - 



(Al) 



where and k locate a given cell in the current block and 
A, B, C, D, x and y are evaluated at cell centres. The x and 
y coordinates are represented by Xij^ and yi,j,h respectively 
(we omit the z coordinate for simplicity). Let the sweep di- 
rection be the x direction, and averaged terms be represented 
by a line over the character. Our averaging technique is as 
followfQ, 

A = 0.5 (Aij,k + Ai-ij^) , (A2) 

where B and C are found in a similar fashion. Derivatives 
in line with the sweep follow: 

d£_ _ ( D i,3,k - A-l,j.fc ' 

dx y Xi^j^k ^i—\,j,h 
Derivatives that are against the sweep follow this form: 
dD _ I (-Dj,j-)-i ; fe — -Di,j-i,fc) + (-D»-i,j+i,fe — Dj-ij-i t k) 

(A4) 

4 Note that we have also used a similar averaging technique which 
averages final terms as a whole, rather than basic values (like 3 X , 
/3ad or in this imaginary case, A). We find no significant difference 
between the two. However, the one presented here more closely 
resembles that which is already employed in the FLASH code. 



(A3) 



F 



A B \ dD dD 
C J dx dy' 



(A5) 



Now we will evaluate the ambipolar diffusion flux terms 
and put them in a form that makes sense for computation. 
In the code we have applied the above averaging technique. 
Please note our matrix notation for flux terms correspond- 
ing to B. We have chosen to make each column correspond 
to B x , By and B z from left to right. Rows for flux and 
source vectors or matrices correspond to the direction of the 
sweep; x direction, y direction and z direction from top to 
bottom. For example, during a sweep in the x direction the 
flux terms only receive additions from the first component 
in the corresponding vector. 

For simplicity we define the vector, 



= -/3 AD B X [B X (V X B)] 
= -Pad 



B X (Byjy + B Z j Z ) 

B y (B x j x + B z j z ) 

B Z (B X j X + Byjy) 



^(Bl+B^) 

jy(B 2 x + B 2 z ) 

j Z {Bl+B 2 y) 



(A6) 



where j = V X B. 

We find that the flux vectors for the magnetic field com- 
ponents are, 



(A7) 




Similarly we define the vector, 
a = (j X B) 

( jyB z -j z B y 
= jzB x - j x B z 

\ fa By — jy B x 



for the source. The sources for ([22} are coded as follows 
terms with a /3ad need to be added, others are already 



(A8) 



(only 
taken 
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care for): 
Sb 



-=-/3ad(V-B) 



a x 











a v 











a z 



(A9) 



where we have split up the source into a vector for each direc- 
tion for simplicity (following our notation outlined above). 
The energy flux works out to, 



would be Az if the sweep was in the z direction. The smallest 
such timestep in a sweep is given to the timestep chooser 
which multiplies the smallest timestep from all processes 
by a CFL factor of typically 0.8. The ambipolar diffusion 
timestep is small and usually dominates over other timesteps 
(such as from the cooling or ideal MHD processes). 



F E += -f3 A oB 2 



(A10) 



The ambipolar diffusion heating term works out to, 

/ al \ 

S E +=/3ad a 2 y , (All) 
V «2 ) 

where again we've chosen how to split up the source. 
Finally, the V ■ B source term works out to, 

/ {B x a x )(V ■ B) \ 
Sb+=-/?ad (B y a y )(V-B) , (A12) 
\ (B z a z )(V-B) J 

where we've chosen to split the source term up partly to save 
coding and match the usage of a's components in Equation 
(|A9I) . The energy equations in this form are quite efficient, 
relying on the calculation of only one component of a for 
each sweep. 

A2 Units, constants and dimensionality 

One advantage of FLASH is that the dimensionless con- 
stants are all 1.0, save the constant for the magnetic field . 
The four principle scaling constants are l|Powell et al.ll 19991 ) : 

rem 

«oo = 1.0 

L s 



1.0 

cm 

L = 1.0 [cm; 



Lcm 3 J 



(A13) 



flQ = 4-7T 



from which most variables can be made dimensionless 
through combinations of these constants. 

Important examples are the magnetic field, B' = 
{a oa ^/Poc,lJ.o)~ 1 B and /3' AD = a 00 p 00 L _1 pio/3AD. Otherwise, 
all other variables are quite straightforward. One can imag- 
ine equations (|20I423(I the same, but drop all jj,o terms to 
1. 



A3 Timesteps 

The timestep is taken from equation (|26|l and written as: 



TAD = Tb- 

= To 



7/AD 



/3ad B 



(A14) 



where To is chosen as ijij m our runs ( as used in 
iMac Low et"ail (|l995T) 1. 

The Ax represents the coordinate of the sweep, so it 
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